markers_for_heatmap <- function(markers) {
res <- NULL
for (i in unique(markers[, 1])) {
tmp <- markers[markers[, 1] == i, ]
if (nrow(tmp) > 10) {
res <- rbind(res, tmp[1:10, ])
} else {
res <- rbind(res, tmp)
}
}
return(res)
}
organise_marker_genes <- function(object, k, p_val, auroc) {
dat <- rowData(object)[, c(paste0("sc3_", k, "_markers_clusts"), paste0("sc3_", k,
"_markers_auroc"), paste0("sc3_", k, "_markers_padj"), "feature_symbol")]
dat <- dat[dat[, paste0("sc3_", k, "_markers_padj")] < p_val & !is.na(dat[, paste0("sc3_",
k, "_markers_padj")]), ]
dat <- dat[dat[, paste0("sc3_", k, "_markers_auroc")] > auroc, ]
d <- NULL
for (i in sort(unique(dat[, paste0("sc3_", k, "_markers_clusts")]))) {
tmp <- dat[dat[, paste0("sc3_", k, "_markers_clusts")] == i, ]
tmp <- tmp[order(tmp[, paste0("sc3_", k, "_markers_auroc")], decreasing = TRUE), ]
d <- rbind(d, tmp)
}
return(d)
}
make_col_ann_for_heatmaps <- function(object, show_pdata) {
if (any(!show_pdata %in% colnames(colData(object)))) {
show_pdata_excl <- show_pdata[!show_pdata %in% colnames(colData(object))]
show_pdata <- show_pdata[show_pdata %in% colnames(colData(object))]
message(paste0("Provided columns '", paste(show_pdata_excl, collapse = "', '"), "' do not exist in the phenoData table!"))
if (length(show_pdata) == 0) {
return(NULL)
}
}
ann <- NULL
if (is.null(metadata(object)$sc3$svm_train_inds)) {
ann <- colData(object)[, colnames(colData(object)) %in% show_pdata]
} else {
ann <- colData(object)[metadata(object)$sc3$svm_train_inds, colnames(colData(object)) %in%
show_pdata]
}
# remove columns with 1 value only
if (length(show_pdata) > 1) {
keep <- unlist(lapply(ann, function(x) {
length(unique(x))
})) > 1
if (!all(keep)) {
message(paste0("Columns '", paste(names(keep)[!keep], collapse = "', '"), "' were excluded from annotation since they contained only a single value."))
}
ann <- ann[, names(keep)[keep]]
if (ncol(ann) == 0) {
ann <- NULL
} else {
ann <- as.data.frame(lapply(ann, function(x) {
if (nlevels(as.factor(x)) > 9)
x else as.factor(x)
}))
# convert outlier scores back to numeric
for (i in grep("_log2_outlier_score", colnames(ann))) {
if (class(ann[, i]) == "factor") {
ann[, i] <- as.numeric(levels(ann[, i]))[ann[, i]]
}
}
}
} else {
if (length(unique(ann)) > 1) {
ann <- as.data.frame(ann)
colnames(ann) <- show_pdata
if (!grepl("_log2_outlier_score", show_pdata)) {
ann <- as.data.frame(lapply(ann, function(x) {
if (nlevels(as.factor(x)) > 9)
return(x) else return(as.factor(x))
}))
}
} else {
message(paste0("Column '", show_pdata, "' was excluded from annotation since they contained only a single value."))
ann <- NULL
}
}
return(ann)
}
sc3_plot_markers.new <- function(object, k, auroc, p.val, show_pdata) {
if (is.null(metadata(object)$sc3$consensus)) {
warning(paste0("Please run sc3_consensus() first!"))
return(object)
}
hc <- metadata(object)$sc3$consensus[[as.character(k)]]$hc
dataset <- get_processed_dataset(object)
if (!is.null(metadata(object)$sc3$svm_train_inds)) {
dataset <- dataset[, metadata(object)$sc3$svm_train_inds]
}
add_ann_col <- FALSE
ann <- NULL
if (!is.null(show_pdata)) {
ann <- make_col_ann_for_heatmaps(object, show_pdata)
if (!is.null(ann)) {
add_ann_col <- TRUE
# make same names for the annotation table
rownames(ann) <- colnames(dataset)
}
}
ann_colors = rep(list(c("grey","white", "firebrick2")),ncol(ann)-sum(sapply(ann, is.factor)))
names(ann_colors) = colnames(ann)[!(sapply(ann, is.factor))]
# get all marker genes
markers <- organise_marker_genes(object, k, p.val, auroc)
# get top 10 marker genes of each cluster
markers <- markers_for_heatmap(markers)
row.ann <- data.frame(Cluster = factor(markers[, 1], levels = unique(markers[, 1])))
rownames(row.ann) <- markers$feature_symbol
do.call(pheatmap::pheatmap, c(list(dataset[markers$feature_symbol, , drop = FALSE], show_colnames = FALSE,
cluster_rows = FALSE, cluster_cols = hc, cutree_cols = k, annotation_row = row.ann, annotation_colors = ann_colors, annotation_names_row = FALSE, cellheight = 10), list(annotation_col = ann)[add_ann_col]))
}
load the data
source("sc_func.R")
suppressMessages(library(Rtsne))
suppressMessages(library(destiny))
suppressMessages(library(ggplot2))
package ‘ggplot2’ was built under R version 3.4.4
suppressMessages(library(gplots))
suppressMessages(library(RColorBrewer))
suppressMessages(library(plotly))
suppressMessages(library(gridExtra))
hmcol<-colorRampPalette(c("blue", "white","red"))(n = 299)
load("sc_192.RData")
do PCA+diffusionmap and extract genes correlated with top PC, excepct PCs that related to cell cycles:
pca_out = prcomp(t(log2(expr_data+1)))
dif <- DiffusionMap( pca_out$x[,c(2,3,5,7,8,9,10,11,12,13,14)],n_eigs = 10)
dif_GO_list = Diff_GO_enrichment(dif,expr_data,p.val=0.001)
Building most specific GOs .....
Loading required package: org.Mm.eg.db
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 2492 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 2097 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 1759 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 1577 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 1093 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 1156 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 1344 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 1594 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 872 nontrivial nodes
parameters:
test statistic: fisher
Building most specific GOs .....
( 9332 GO terms found. )
Build GO DAG topology ..........
( 13302 GO terms and 31005 relations. )
Annotating nodes ...............
( 9474 genes annotated to the GO terms. )
-- Classic Algorithm --
the algorithm is scoring 762 nontrivial nodes
parameters:
test statistic: fisher
diff_genes = c()
for(i in 1:10){diff_genes = c(diff_genes,rownames(dif_GO_list$DE_genes[[i]]))}
DE_expr = expr_data[rownames(expr_data) %in% diff_genes,]
DE_clean = DE_expr[apply(log2(DE_expr+1),1,mean)<6,]
DE_clean = DE_clean[rowMeans(DE_clean)>0.6,]
dim(DE_expr)
[1] 407 187
hr <- hclust(dist(cor(log2(DE_clean+1))))
tsne_ph = Rtsne(log10(ph_data[,!(colnames(ph_data) %in% c("CD11b"))]),pca = FALSE,perplexity = 10)
tsne_ph1 = Rtsne(log10(ph_data[,!(colnames(ph_data) %in% c("CD11b"))]),dim=1,pca = FALSE,perplexity = 10)
# define some clusters
mycl <- cutree(hr, k=10)
table(mycl)
mycl
1 2 3 4 5 6 7 8 9 10
10 4 42 31 46 35 6 2 10 1
#mycl[mycl>8] = 8 # exclude rare cluster
# get a color palette equal to the number of clusters
#clusterCols <- rainbow(length(unique(mycl)))
clusterCols = brewer.pal(10, "Set3")
# create vector of colors for side bar
myClusterSideBar <- clusterCols[mycl]
#ph_data$cluster = myClusterSideBar
DE_hm = heatmap.2(log2(DE_clean[,order(mycl)]+1),
Colv=FALSE,
trace="none",
dendrogram = "none",
col=hmcol,
density.info = "none",
labCol = "",
labRow = "",
ColSideColors = myClusterSideBar[order(mycl)],
key = FALSE
)
ma = log2(DE_clean[,order(mycl)]+1)
ma = ma[DE_hm$rowInd,]
plot_ly(
x = colnames(ma), y = rownames(ma),
z = ma, type = "heatmap",colorscale = hmcol
)
database = "immgen"
immgen_samples <- read.delim(paste0("/Users/tian.l/Dropbox/research/Dach1_paper/ref_data/",database,"_samples.txt"), stringsAsFactors=FALSE)
immgen_probes <- read.delim(paste0("/Users/tian.l/Dropbox/research/Dach1_paper/ref_data/",database,"_probes.txt"), header=FALSE, stringsAsFactors=FALSE)
immgen_expression <- read.delim(paste0("/Users/tian.l/Dropbox/research/Dach1_paper/ref_data/",database,"_expression.txt"), stringsAsFactors=FALSE)
library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id", "external_gene_name", "description"),values=immgen_probes$V2 ,mart=mart)
Batch submitting query [====-------------------------------------------------------------------------------] 4% eta: 29s
Batch submitting query [======-----------------------------------------------------------------------------] 7% eta: 29s
Batch submitting query [=======----------------------------------------------------------------------------] 9% eta: 26s
Batch submitting query [=========--------------------------------------------------------------------------] 11% eta: 24s
Batch submitting query [===========------------------------------------------------------------------------] 13% eta: 22s
Batch submitting query [=============----------------------------------------------------------------------] 16% eta: 21s
Batch submitting query [===============--------------------------------------------------------------------] 18% eta: 20s
Batch submitting query [=================------------------------------------------------------------------] 20% eta: 19s
Batch submitting query [==================-----------------------------------------------------------------] 22% eta: 18s
Batch submitting query [====================---------------------------------------------------------------] 24% eta: 17s
Batch submitting query [======================-------------------------------------------------------------] 27% eta: 17s
Batch submitting query [========================-----------------------------------------------------------] 29% eta: 17s
Batch submitting query [==========================---------------------------------------------------------] 31% eta: 17s
Batch submitting query [============================-------------------------------------------------------] 33% eta: 17s
Batch submitting query [==============================-----------------------------------------------------] 36% eta: 17s
Batch submitting query [===============================----------------------------------------------------] 38% eta: 16s
Batch submitting query [=================================--------------------------------------------------] 40% eta: 15s
Batch submitting query [===================================------------------------------------------------] 42% eta: 14s
Batch submitting query [=====================================----------------------------------------------] 44% eta: 14s
Batch submitting query [=======================================--------------------------------------------] 47% eta: 13s
Batch submitting query [=========================================------------------------------------------] 49% eta: 12s
Batch submitting query [==========================================-----------------------------------------] 51% eta: 12s
Batch submitting query [============================================---------------------------------------] 53% eta: 11s
Batch submitting query [==============================================-------------------------------------] 56% eta: 11s
Batch submitting query [================================================-----------------------------------] 58% eta: 10s
Batch submitting query [==================================================---------------------------------] 60% eta: 10s
Batch submitting query [====================================================-------------------------------] 62% eta: 9s
Batch submitting query [=====================================================------------------------------] 64% eta: 8s
Batch submitting query [=======================================================----------------------------] 67% eta: 8s
Batch submitting query [=========================================================--------------------------] 69% eta: 7s
Batch submitting query [===========================================================------------------------] 71% eta: 7s
Batch submitting query [=============================================================----------------------] 73% eta: 6s
Batch submitting query [===============================================================--------------------] 76% eta: 6s
Batch submitting query [=================================================================------------------] 78% eta: 5s
Batch submitting query [==================================================================-----------------] 80% eta: 5s
Batch submitting query [====================================================================---------------] 82% eta: 4s
Batch submitting query [======================================================================-------------] 84% eta: 4s
Batch submitting query [========================================================================-----------] 87% eta: 3s
Batch submitting query [==========================================================================---------] 89% eta: 3s
Batch submitting query [============================================================================-------] 91% eta: 2s
Batch submitting query [=============================================================================------] 93% eta: 2s
Batch submitting query [===============================================================================----] 96% eta: 1s
Batch submitting query [=================================================================================--] 98% eta: 1s
Batch submitting query [===================================================================================] 100% eta: 0s
immgen_probes = immgen_probes[immgen_probes$V2 %in% G_list$ensembl_gene_id,]
G_list = G_list[match(immgen_probes$V2, G_list$ensembl_gene_id),]
immgen_probes$external_gene_name = G_list$external_gene_name
immgen_expression = immgen_expression[immgen_expression$probeId %in% immgen_probes$V1,]
immgen_probes = immgen_probes[immgen_probes$V1 %in% immgen_expression$probeId,]
immgen_probes = immgen_probes[match(immgen_expression$probeId, immgen_probes$V1),]
immgen_expression$gene_id = immgen_probes$external_gene_name
immgen_expression = immgen_expression[!duplicated(immgen_expression$gene_id),]
g_i = immgen_expression$gene_id
immgen_expression = immgen_expression[,!(colnames(immgen_expression) %in% c("probeId","gene_id"))]
immgen_expression = as.matrix(immgen_expression)
rownames(immgen_expression) = g_i
table(immgen_samples$cell_lineage)
ab T cells B cells Dendritic cells gd T cells Innate Lymphocytes Macrophages
67 26 36 22 22 15
Monocytes Neutrophils Stem Cells Stromal cells
8 6 14 8
immgen_samples$sampleId = gsub("-",".",immgen_samples$sampleId)
#SC_list = immgen_samples[immgen_samples$cell_lineage == "Stem Cells", "sampleId"]
#SC_list = immgen_samples$sampleId
SC_list = c("SC_CDP_BM","SC_CMP_BM","SC_GMP_BM",
"SC_LT34F_BM","SC_LTSL_BM","SC_MDP_BM",
"SC_MEP_BM","SC_MPP34F_BM","SC_ST34F_BM",
"SC_STSL_BM","proB_CLP_BM","preT_ETP_Th")
sub_immgen_expression = immgen_expression[rownames(immgen_expression) %in% rownames(DE_expr), colnames(immgen_expression) %in% SC_list]
sub_immgen_expression = sub_immgen_expression[!duplicated(rownames(sub_immgen_expression)),]
sub_DE_expr = DE_expr[rownames(DE_expr) %in% rownames(sub_immgen_expression),]
sub_DE_expr = sub_DE_expr[!duplicated(rownames(sub_DE_expr)),]
sub_DE_expr = sub_DE_expr[order(rownames(sub_DE_expr)),]
sub_immgen_expression = sub_immgen_expression[order(rownames(sub_immgen_expression)),]
p.val_mat = c()
for (i in 1:ncol(sub_immgen_expression))
{
p.val_vec = apply(sub_DE_expr, 2, function(x){
cor.test(x,sub_immgen_expression[,i],
method = "spearman",
alternative = "greater",exact=FALSE)$p.value})
p.val_mat = rbind(p.val_mat,-log10(p.val_vec))
}
rownames(p.val_mat) = colnames(sub_immgen_expression)
hm_immgen = heatmap.2(p.val_mat[,order(mycl)],trace="none",
dendrogram="none",
Colv=FALSE,
col=hmcol,
scale="column",
labRow = "",
#ColSideColors = myClusterSideBar[order(mycl)],
labCol = "",
density.info="none",
key=FALSE)
p_ma = p.val_mat[,order(mycl)]
p_ma = p_ma[hm_immgen$rowInd,]
m = list(
l = 100,
r = 40,
b = 10,
t = 10,
pad = 0
)
plot_ly(
x = colnames(p_ma), y = rownames(p_ma),
z = p_ma, type = "heatmap",colorscale = hmcol
)%>%
layout(autosize = F, margin = m)
combine two heatmap:
library(ComplexHeatmap)
Loading required package: grid
Attaching package: ‘grid’
The following object is masked from ‘package:topGO’:
depth
========================================
ComplexHeatmap version 1.17.1
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://bioconductor.org/packages/ComplexHeatmap/
If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
genomic data. Bioinformatics 2016.
========================================
Attaching package: ‘ComplexHeatmap’
The following object is masked from ‘package:plotly’:
add_heatmap
ht1 = Heatmap(t(p.val_mat), name = "Immgen", show_row_names=FALSE, width = 1)
ht2 = Heatmap(t(log2(DE_clean+1)), name = "RNA-seq", show_row_names=FALSE, show_column_names=FALSE, width = 2)
ht3 = Heatmap(scale(log10(ph_data[,!(colnames(ph_data)=="CD11b")]+1)), name = "FACS",show_row_names=FALSE, width = 0.5)
ht_list = draw(ht3+ht2+ht1, main_heatmap = "RNA-seq", km = 10,show_row_dend=FALSE)
#clusterCols = brewer.pal(10, "Set3")
# create vector of colors for side bar
#myClusterSideBar <- clusterCols[mycl]
ht_list
tsne_diffmap_DE = Rtsne(cor(log2(DE_expr+1)))
tsne_diffmap_DE1 = Rtsne(cor(log2(DE_expr+1)),dim=1)
plot(tsne_diffmap_DE$Y[,1], tsne_diffmap_DE$Y[,2])
plot(tsne_diffmap_DE1$Y, tsne_ph1$Y)
cor_shalin =cor(t(expr_data),method="spearman")
cor_Dach1 = cor_shalin["Dach1",]
cor_Dach1 = cor_Dach1[order(cor_Dach1)]
genes that positively correlated with Dach1:
as.data.frame(head(cor_Dach1[order(cor_Dach1,decreasing=TRUE)],n=50))
some scatter plot:
g_e = as.data.frame(t(log2(expr_data+1)))
p1=ggplot(data=g_e,aes(y=Dach1, x=Muc13))+geom_point()
p2=ggplot(data=g_e,aes(y=Dach1, x=Gata2))+geom_point()
p3=ggplot(data=g_e,aes(y=Dach1, x=Cdk17))+geom_point()
p4=ggplot(data=g_e,aes(y=Dach1, x=Sox6))+geom_point()
grid.arrange(p1, p2, p3, p4, ncol=2)
Error: `data` must be uniquely named but has duplicate elements
genes that negatively correlated with Dach1:
some scatter plot:
p1=ggplot(data=g_e,aes(y=Dach1, x=Evl))+geom_point()
p2=ggplot(data=g_e,aes(y=Dach1, x=Il12a))+geom_point()
p3=ggplot(data=g_e,aes(y=Dach1, x=Dntt))+geom_point()
p4=ggplot(data=g_e,aes(y=Dach1, x=Flt3))+geom_point()
grid.arrange(p1, p2, p3, p4, ncol=2)
Error: `data` must be uniquely named but has duplicate elements
data from Gottgens’s Blood paper, we only select cells labeled as “HSC” and “LT-HSC”.
coordinates_gene_counts_flow_cytometry <- read.delim("/Users/tian.l/Dropbox/research/Dach1_paper/public_data/coordinates_gene_counts_flow_cytometry.txt", row.names=1, stringsAsFactors=FALSE)
coordinates_gene_counts_flow_cytometry = coordinates_gene_counts_flow_cytometry[complete.cases(coordinates_gene_counts_flow_cytometry),]
gene_expr = coordinates_gene_counts_flow_cytometry[,-(1:26)]
phen_data = coordinates_gene_counts_flow_cytometry[,6:15]
gene_expr_t = t(gene_expr)
G_cor = cor(t(gene_expr_t))
library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id", "external_gene_name", "description"),values=rownames(gene_expr_t) ,mart= mart)
gene_expr_t = gene_expr_t[rownames(gene_expr_t) %in% G_list$ensembl_gene_id,]
G_list = G_list[match(rownames(gene_expr_t), G_list$ensembl_gene_id),]
rownames(gene_expr_t) = G_list$external_gene_name
gene_expr_sub = gene_expr_t[,coordinates_gene_counts_flow_cytometry$group %in% c("HSC", "LT-HSC")]
corr_sub = cor(t(gene_expr_sub), method="spearman")
cor_Dach1_pub = corr_sub["Dach1",]
cor_Dach1_pub = cor_Dach1_pub[order(cor_Dach1_pub)]
genes that positively correlated with Dach1:
as.data.frame(head(cor_Dach1_pub[order(cor_Dach1_pub,decreasing=TRUE)],n=50))
genes that negatively correlated with Dach1:
as.data.frame(head(cor_Dach1_pub[order(cor_Dach1_pub)],n=50))
Genes that negatively correlated to Dach1 in both dataset (first 50):
neg = toupper(names(head(cor_Dach1_pub[order(cor_Dach1_pub)],n=50)))
neg = neg[neg %in% toupper(names(head(cor_Dach1[order(cor_Dach1)],n=50)))]
neg
Genes that positively correlated to Dach1 in both dataset (first 50):
pos = toupper(names(head(cor_Dach1_pub[order(cor_Dach1_pub,decreasing=TRUE)],n=50)))
pos = pos[pos %in% toupper(names(head(cor_Dach1[order(cor_Dach1,decreasing=TRUE)],n=50)))]
pos
Cells are classified as being in G1 phase if the G1 score is above 0.5 and greater than the G2/M score; in G2/M phase if the G2/M score is above 0.5 and greater than the G1 score; and in S phase if neither score is above 0.5.
library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
library(org.Mm.eg.db)
anno <- select(org.Mm.eg.db, keys=rownames(expr_data), keytype="SYMBOL", column="ENSEMBL")
ensembl <- anno$ENSEMBL[match(rownames(expr_data), anno$SYMBOL)]
assignments <- cyclone(log2(expr_data+1), mm.pairs, gene.names=ensembl)
ggplot(data=NULL, aes(x=assignments$score$G1, y=assignments$score$G2M, col=log2(expr_data["Dach1",]+1)))+
geom_point()+
geom_vline(xintercept = .5)+
geom_hline(yintercept = .5)
cell_state = rep("S",ncol(expr_data))
cell_state[assignments$score$G1>0.5] = "G0/G1"
cell_state[assignments$score$G2M>0.5] = "G2/M"
dach1 = log2(expr_data["Dach1",]+1)
ggplot(data=NULL,aes(x=cell_state,y=dach1,fill=cell_state))+
geom_violin(alpha=0.7) + geom_jitter(height = 0, width = 0.1)+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
public_assignments <- cyclone(t(gene_expr)[,coordinates_gene_counts_flow_cytometry$group %in% c("HSC", "LT-HSC")], mm.pairs, gene.names=rownames(t(gene_expr)))
ggplot(data=NULL, aes(x=public_assignments$score$G1, y=public_assignments$score$G2M, col=gene_expr_t["Dach1",coordinates_gene_counts_flow_cytometry$group %in% c("HSC", "LT-HSC")]))+
geom_point()+
labs(x="G1",y="G2M",color="Dach1 expr")+
geom_vline(xintercept = .5)+
geom_hline(yintercept = .5)
cell_state = rep("S",sum(coordinates_gene_counts_flow_cytometry$group %in% c("HSC", "LT-HSC")))
cell_state[public_assignments$score$G1>0.5] = "G0/G1"
cell_state[public_assignments$score$G2M>0.5] = "G2/M"
dach1 = gene_expr_t["Dach1",]
ggplot(data=NULL,aes(x=cell_state,y=dach1,fill=cell_state))+
geom_violin(alpha=0.7) + geom_jitter(height = 0, width = 0.1)+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
sceset <- sc3(sceset, ks = 4:8, biology = TRUE,pct_dropout_max=95,rand_seed=19910603,kmeans_nstart=50)
Setting SC3 parameters...
Calculating distances between the cells...
starting worker pid=20265 on localhost:11927 at 11:29:35.450
starting worker pid=20273 on localhost:11927 at 11:29:35.660
starting worker pid=20281 on localhost:11927 at 11:29:35.868
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
loaded SC3 and set parent environment
loaded SC3 and set parent environment
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: foreach
Loading required package: foreach
Loading required package: rngtools
Loading required package: rngtools
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: pkgmaker
Loading required package: pkgmaker
Loading required package: registry
Loading required package: registry
Loading required package: registry
Attaching package: ‘pkgmaker’
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
The following object is masked from ‘package:base’:
isNamespaceLoaded
Performing transformations and calculating eigenvectors...
starting worker pid=20293 on localhost:11927 at 11:29:41.435
starting worker pid=20301 on localhost:11927 at 11:29:41.666
starting worker pid=20309 on localhost:11927 at 11:29:41.893
starting worker pid=20317 on localhost:11927 at 11:29:42.121
starting worker pid=20325 on localhost:11927 at 11:29:42.349
starting worker pid=20333 on localhost:11927 at 11:29:42.574
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: rngtools
Loading required package: registry
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Performing k-means clustering...
starting worker pid=20348 on localhost:11927 at 11:29:50.050
starting worker pid=20356 on localhost:11927 at 11:29:50.276
starting worker pid=20364 on localhost:11927 at 11:29:50.501
starting worker pid=20372 on localhost:11927 at 11:29:50.729
starting worker pid=20380 on localhost:11927 at 11:29:50.956
starting worker pid=20388 on localhost:11927 at 11:29:51.191
starting worker pid=20396 on localhost:11927 at 11:29:51.417
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
loaded SC3 and set parent environment
Loading required package: pkgmaker
Loading required package: foreach
Loading required package: registry
loaded SC3 and set parent environment
Loading required package: rngtools
Loading required package: foreach
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Loading required package: pkgmaker
Loading required package: registry
Loading required package: rngtools
Loading required package: pkgmaker
loaded SC3 and set parent environment
Attaching package: ‘pkgmaker’
Loading required package: registry
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: foreach
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Loading required package: rngtools
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: pkgmaker
Loading required package: registry
Loading required package: registry
loaded SC3 and set parent environment
Loading required package: foreach
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
|
|
|= | 0%
| | 0%
|
|= | 1%
|
|== | 1%
|
|
|=== | 3%
|
|== | 2%
|=== | 2%
|
|==== | 3%
|
|==== | 3%
|
|===== | 4%
|
|===== | 4%
|
|====== | 5%
|
|====== | 5%
|
|======= | 5%
|
|======= | 6%
|
|======== | 6%
|
|======== | 7%
|
|========= | 7%
|
|========= | 8%
|
|========== | 8%
|
|========== | 8%
|
|=========== | 9%
|
|=========== | 9%
|
|============ | 10%
|
|============ | 10%
|
|============= | 10%
|
|============= | 11%
|
|============== | 11%
|
|============== | 12%
|
|=============== | 12%
|
|=============== | 13%
|
|================ | 13%
|
|================ | 13%
|
|================= | 14%
|
|================= | 14%
|
|================== | 15%
|
|================== | 15%
|
|=================== | 15%
|
|=================== | 16%
|
|==================== | 16%
|
|==================== | 17%
|
|===================== | 17%
|
|===================== | 18%
|
|====================== | 18%
|
|====================== | 18%
|
|======================= | 19%
|
|======================= | 19%
|
|======================== | 20%
|
|========================= | 20%
|
|========================= | 21%
|
|========================== | 21%
|
|========================== | 21%
|
|=========================== | 22%
|
|=========================== | 22%
|
|============================ | 23%
|
|============================ | 23%
|
|============================= | 23%
|
|============================= | 24%
|
|============================== | 24%
|
|============================== | 25%
|
|=============================== | 25%
|
|=============================== | 26%
|
|================================ | 26%
|
|================================ | 26%
|
|================================= | 27%
|
|================================= | 27%
|
|================================== | 28%
|
|================================== | 28%
|
|=================================== | 28%
|
|=================================== | 29%
|
|==================================== | 29%
|
|==================================== | 30%
|
|===================================== | 30%
|
|===================================== | 31%
|
|====================================== | 31%
|
|====================================== | 31%
|
|======================================= | 32%
|
|======================================= | 32%
|
|======================================== | 33%
|
|======================================== | 33%
|
|========================================= | 33%
|
|========================================= | 34%
|
|========================================== | 34%
|
|========================================== | 35%
|
|=========================================== | 35%
|
|=========================================== | 36%
|
|============================================ | 36%
|
|============================================ | 36%
|
|============================================= | 37%
|
|============================================= | 37%
|
|============================================== | 38%
|
|============================================== | 38%
|
|=============================================== | 38%
|
|=============================================== | 39%
|
|================================================ | 39%
|
|================================================ | 40%
|
|================================================= | 40%
|
|================================================== | 41%
|
|================================================== | 41%
|
|=================================================== | 41%
|
|=================================================== | 42%
|
|==================================================== | 42%
|
|==================================================== | 43%
|
|===================================================== | 43%
|
|===================================================== | 44%
|
|====================================================== | 44%
|
|====================================================== | 44%
|
|======================================================= | 45%
|
|======================================================= | 45%
|
|======================================================== | 46%
|
|======================================================== | 46%
|
|========================================================= | 46%
|
|========================================================= | 47%
|
|========================================================== | 47%
|
|========================================================== | 48%
|
|=========================================================== | 48%
|
|=========================================================== | 49%
|
|============================================================ | 49%
|
|============================================================ | 49%
|
|============================================================= | 50%
|
|============================================================= | 50%
|
|============================================================== | 51%
|
|============================================================== | 51%
|
|=============================================================== | 51%
|
|=============================================================== | 52%
|
|================================================================ | 52%
|
|================================================================ | 53%
|
|================================================================= | 53%
|
|================================================================= | 54%
|
|================================================================== | 54%
|
|================================================================== | 54%
|
|=================================================================== | 55%
|
|=================================================================== | 55%
|
|==================================================================== | 56%
|
|==================================================================== | 56%
|
|===================================================================== | 56%
|
|===================================================================== | 57%
|
|====================================================================== | 57%
|
|====================================================================== | 58%
|
|======================================================================= | 58%
|
|======================================================================= | 59%
|
|======================================================================== | 59%
|
|======================================================================== | 59%
|
|========================================================================= | 60%
|
|========================================================================== | 60%
|
|========================================================================== | 61%
|
|=========================================================================== | 61%
|
|=========================================================================== | 62%
|
|============================================================================ | 62%
|
|============================================================================ | 62%
|
|============================================================================= | 63%
|
|============================================================================= | 63%
|
|============================================================================== | 64%
|
|============================================================================== | 64%
|
|=============================================================================== | 64%
|
|=============================================================================== | 65%
|
|================================================================================ | 65%
|
|================================================================================ | 66%
|
|================================================================================= | 66%
|
|================================================================================= | 67%
|
|================================================================================== | 67%
|
|================================================================================== | 67%
|
|=================================================================================== | 68%
|
|=================================================================================== | 68%
|
|==================================================================================== | 69%
|
|==================================================================================== | 69%
|
|===================================================================================== | 69%
|
|===================================================================================== | 70%
|
|====================================================================================== | 70%
|
|====================================================================================== | 71%
|
|======================================================================================= | 71%
|
|======================================================================================= | 72%
|
|======================================================================================== | 72%
|
|======================================================================================== | 72%
|
|========================================================================================= | 73%
|
|========================================================================================= | 73%
|
|========================================================================================== | 74%
|
|========================================================================================== | 74%
|
|=========================================================================================== | 74%
|
|=========================================================================================== | 75%
|
|============================================================================================ | 75%
|
|============================================================================================ | 76%
|
|============================================================================================= | 76%
|
|============================================================================================= | 77%
|
|============================================================================================== | 77%
|
|============================================================================================== | 77%
|
|=============================================================================================== | 78%
|
|=============================================================================================== | 78%
|
|================================================================================================ | 79%
|
|================================================================================================ | 79%
|
|================================================================================================= | 79%
|
|================================================================================================= | 80%
|
|================================================================================================== | 80%
|
|=================================================================================================== | 81%
|
|=================================================================================================== | 81%
|
|==================================================================================================== | 82%
|
|==================================================================================================== | 82%
|
|===================================================================================================== | 82%
|
|===================================================================================================== | 83%
|
|====================================================================================================== | 83%
|
|====================================================================================================== | 84%
|
|======================================================================================================= | 84%
|
|======================================================================================================= | 85%
|
|======================================================================================================== | 85%
|
|======================================================================================================== | 85%
|
|========================================================================================================= | 86%
|
|========================================================================================================= | 86%
|
|========================================================================================================== | 87%
|
|========================================================================================================== | 87%
|
|=========================================================================================================== | 87%
|
|=========================================================================================================== | 88%
|
|============================================================================================================ | 88%
|
|============================================================================================================ | 89%
|
|============================================================================================================= | 89%
|
|============================================================================================================= | 90%
|
|============================================================================================================== | 90%
|
|============================================================================================================== | 90%
|
|=============================================================================================================== | 91%
|
|=============================================================================================================== | 91%
|
|================================================================================================================ | 92%
|
|================================================================================================================ | 92%
|
|================================================================================================================= | 92%
|
|================================================================================================================= | 93%
|
|================================================================================================================== | 93%
|
|================================================================================================================== | 94%
|
|=================================================================================================================== | 94%
|
|=================================================================================================================== | 95%
|
|==================================================================================================================== | 95%
|
|==================================================================================================================== | 95%
|
|===================================================================================================================== | 96%
|
|===================================================================================================================== | 96%
|
|====================================================================================================================== | 97%
|
|====================================================================================================================== | 97%
|
|======================================================================================================================= | 97%
|
|======================================================================================================================= | 98%
|
|======================================================================================================================== | 98%
|
|======================================================================================================================== | 99%
|
|========================================================================================================================= | 99%
|
|========================================================================================================================= | 100%
|
|==========================================================================================================================| 100%
Calculating consensus matrix...
starting worker pid=20413 on localhost:11927 at 11:30:01.318
starting worker pid=20421 on localhost:11927 at 11:30:01.554
starting worker pid=20429 on localhost:11927 at 11:30:01.791
starting worker pid=20437 on localhost:11927 at 11:30:02.021
starting worker pid=20445 on localhost:11927 at 11:30:02.248
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
loaded SC3 and set parent environment
Attaching package: ‘pkgmaker’
Loading required package: foreach
The following object is masked from ‘package:base’:
isNamespaceLoaded
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
loaded SC3 and set parent environment
loaded SC3 and set parent environment
Loading required package: registry
Loading required package: foreach
Loading required package: foreach
Loading required package: rngtools
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: pkgmaker
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Loading required package: registry
Loading required package: registry
Attaching package: ‘pkgmaker’
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
The following object is masked from ‘package:base’:
isNamespaceLoaded
Calculating biology...
starting worker pid=20459 on localhost:11927 at 11:30:08.867
starting worker pid=20469 on localhost:11927 at 11:30:09.107
starting worker pid=20477 on localhost:11927 at 11:30:09.333
starting worker pid=20485 on localhost:11927 at 11:30:09.560
starting worker pid=20493 on localhost:11927 at 11:30:09.789
starting worker pid=20501 on localhost:11927 at 11:30:10.021
starting worker pid=20509 on localhost:11927 at 11:30:10.254
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
Loading required package: SC3
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
loaded SC3 and set parent environment
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: pkgmaker
Loading required package: foreach
Loading required package: registry
Loading required package: rngtools
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: pkgmaker
Loading required package: rngtools
loaded SC3 and set parent environment
Loading required package: registry
Loading required package: registry
Loading required package: pkgmaker
Loading required package: foreach
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Loading required package: rngtools
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
loaded SC3 and set parent environment
Loading required package: foreach
Loading required package: rngtools
Loading required package: pkgmaker
Loading required package: registry
Attaching package: ‘pkgmaker’
The following object is masked from ‘package:base’:
isNamespaceLoaded
No outliers detected in cluster 1. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 1. Distribution of gene expression in cells is too skewed towards 0.
The covariance matrix has become singular during
the iterations of the MCD algorithm.
There are 263 observations (in the entire dataset of 323 obs.) lying on
the hyperplane with equation a_1*(x_i1 - m_1) + ... + a_p*(x_ip - m_p)
= 0 with (m_1, ..., m_p) the mean of these observations and
coefficients a_i from the vector a <- c(0.3565493, 0.5424786,
0.3488608, -0.5984084, -0.3079517, 0.0590775, -0.0216488)The covariance matrix has become singular during
the iterations of the MCD algorithm.
There are 0 observations (in the entire dataset of 323 obs.) lying on
the hyperplane with equation a_1*(x_i1 - m_1) + ... + a_p*(x_ip - m_p)
= 0 with (m_1, ..., m_p) the mean of these observations and
coefficients a_i from the vector a <- c(-0.3254948, 0.5212834,
-0.6797406, 0.0837222, -0.3883333, 0.0299857, 0.0394742)No outliers detected in cluster 3. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 2. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 2. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 2. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 3. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 3. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 4. Distribution of gene expression in cells is too skewed towards 0.
The covariance matrix has become singular during
the iterations of the MCD algorithm.
There are 0 observations (in the entire dataset of 323 obs.) lying on
the hyperplane with equation a_1*(x_i1 - m_1) + ... + a_p*(x_ip - m_p)
= 0 with (m_1, ..., m_p) the mean of these observations and
coefficients a_i from the vector a <- c(0.3254948, -0.5212834,
0.6797406, -0.0837222, 0.3883333, -0.0299857, -0.0394742)No outliers detected in cluster 6. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 7. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 5. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 5. Distribution of gene expression in cells is too skewed towards 0.
No outliers detected in cluster 8. Distribution of gene expression in cells is too skewed towards 0.
sceset <- sc3_estimate_k(sceset)
Estimating k...
sceset$SC_LTSL_BM = p.val_mat["SC_LTSL_BM",]
sceset$SC_STSL_BM = p.val_mat["SC_STSL_BM",]
sceset$SC_MEP_BM = p.val_mat["SC_MEP_BM",]
sceset$SC_CMP_BM = p.val_mat["SC_CMP_BM",]
sceset$proB_CLP_BM = p.val_mat["proB_CLP_BM",]
sceset$SC_GMP_BM = p.val_mat["SC_GMP_BM",]
sceset$CD150 = log2(ph_data[, "CD150"]+1)
sceset$CD135 = log2(ph_data[, "CD135"]+1)
sceset$CD135[(sceset)$CD135<5] = 5
sceset$CD127 = log2(ph_data[, "CD127"]+1)
sceset$CD127[(sceset)$CD127<5] = 5
sceset$CD16.32 = log2(ph_data[, "CD16.32"]+1)
sceset$Sca1 = log2(ph_data[, "Sca.1"]+1)
sceset$cKit = log2(ph_data[, "cKit"]+1)
#sceset$G1 = assignments$scores$G1
#sceset$G2M = assignments$scores$G2M
pdf("Dach1_clusters_IMMGEN.pdf",width = 8,height = 10)
sc3_plot_markers(sceset, k = 7,show_pdata = c(
"SC_LTSL_BM",
"SC_STSL_BM",
"SC_MEP_BM",
"SC_CMP_BM",
"proB_CLP_BM",
"SC_GMP_BM",
"sc3_7_clusters"
), p.val=0.01,auroc=0.75)
dev.off()
pdf("Dach1_clusters_FACs.pdf",width = 8,height = 9)
sc3_plot_markers(sceset, k = 7,show_pdata = c(
"CD150",
"CD135",
"CD127",
"CD16.32",
"G1",
"G2M"
), p.val=0.01,auroc=0.75)
dev.off()
pdf("Dach1_clusters_all.pdf",width = 8,height = 16)
sc3_plot_markers.new(sceset, k = 5,show_pdata = c(
"CD150",
"CD135",
"CD127",
"CD16.32",
"SC_LTSL_BM",
"SC_STSL_BM",
"SC_MEP_BM",
"SC_CMP_BM",
"proB_CLP_BM",
"SC_GMP_BM"
), p.val=0.05,auroc=0.6)
dev.off()
RStudioGD
2
library(edgeR)
y = DGEList(counts=expr_data)
y = estimateCommonDisp(y)
design <- model.matrix(~log(expr_data["Dach1",]+1))
fit = glmFit(y,design=design)
lrt = glmLRT(fit,coef = 2)
t = topTags(lrt,n=Inf)
test_df = t@.Data[[1]]
test_df = test_df[!(rownames(test_df)=="Dach1"),]
test_df$FDR = -log10(test_df$FDR)
plot(test_df$logFC, test_df$FDR)
test_df["Gata2",]
library(scran)
sceall <- newSCESet(countData=expr_data[!duplicated(rownames(expr_data)),])
sceall = calculateQCMetrics(sceall)
var.fit <- trendVar(sceall, trend="loess", use.spikes=FALSE, span=0.2)
var.out <- decomposeVar(sceall, var.fit)
hvg.out <- var.out[order(var.out$bio, decreasing=TRUE)[1:2000],]
scran_cor = correlatePairs(sceall, subset.row=rownames(hvg.out),BPPARAM=MulticoreParam(workers = 4))
save(scran_cor,file="/Users/tian.l/Dropbox/research/Dach1_paper/scran_cor_top2000.RData")
scran_cor_dach1 = scran_cor[scran_cor$gene1 == "Dach1" | scran_cor$gene2 == "Dach1",]
the_cor_genes = rep(NA, nrow(scran_cor_dach1))
the_cor_genes[scran_cor_dach1$gene1 == "Dach1"] = scran_cor_dach1$gene2[scran_cor_dach1$gene1 == "Dach1"]
the_cor_genes[scran_cor_dach1$gene2 == "Dach1"] = scran_cor_dach1$gene1[scran_cor_dach1$gene2 == "Dach1"]
library(readr)
renameME <- read_csv("~/Dropbox/research/Dach1_paper/renameME_LTMPP.csv")
renameME_DE = renameME[renameME$adj.P.Val<0.05,]
renameME_DE_CMP = renameME_DE[renameME_DE$FC>1,]
renameME_DE_CLP = renameME_DE[renameME_DE$FC<1,]
tmp = scran_cor_dach1$rho[scran_cor_dach1$FDR<0.05]
pdf("barcode_plot_Dach1_immgen_LT_MPP.pdf")
barcodeplot(scran_cor_dach1$rho, index=(the_cor_genes %in% renameME_DE_CMP$Gene_Symbol ),
index2=(the_cor_genes %in% renameME_DE_CLP$Gene_Symbol ), main="LTvsMPP34",
labels=c("negative correlation","positive correlation"),
quantiles=c(max(tmp[tmp<0]),min(tmp[tmp>0])),
xlab="speaman correlation to Dach1 expressions"
)
dev.off()
hm20 = colorRampPalette(c("blue", "white","red"))(n = 20)
side_bar = c("orchid1", "steelblue1")[as.factor(scran_cor_dach1$rho[scran_cor_dach1$FDR<0.05]>0)]
names(side_bar) = "rho to Dach1"
top_dach1_DE = the_cor_genes[scran_cor_dach1$FDR<0.05]
pdf("top_correlated_to_Dach1.pdf")
heatmap.2(cor(t(log2(expr_data[top_dach1_DE,]+1)),method = "spearman"),
trace = "none",
col=hmcol,
main="Correlation of genes correlated with Dach1",
ColSideColors=side_bar,
RowSideColors=side_bar,
key.title="speaman rho")
legend("topright",
legend = c("positive correlated to Dach1", "negatively correlated to Dach1"),
col = c("steelblue1","orchid1"),
lty= 1.5, lwd = 2,
cex=.6)
dev.off()
CODEX[CODEX$Factor %in% the_cor_genes,]
scran_cor_dach1[the_cor_genes %in% CODEX$Factor,]
TF_genes = scPipe::get_genes_by_GO(returns="external_gene_name",
dataset="mmusculus_gene_ensembl",
go="GO:0003700")
pdf("FACS_overlay.pdf",width = 4,height = 3.5)
# cluster = as.factor(sceset$sc3_5_clusters == 1)
# ggplot(data=as.data.frame(colData(sceset)),aes(x=CD135,y=CD150,col=cluster))+
# theme_bw()+
# geom_point()+
# scale_color_manual(values=c("grey80","red"))+
# #labs(color="stem")+
# ggtitle("stem")
#
# cluster = as.factor(sceset$sc3_5_clusters %in% c(1,2))
# ggplot(data=as.data.frame(colData(sceset)),aes(x=CD135,y=CD150,col=cluster))+
# theme_bw()+
# geom_point()+
# scale_color_manual(values=c("grey80","red"))+
# #labs(color="stem+MEP")+
# ggtitle("stem+MEP")
#
# cluster = as.factor(sceset$sc3_5_clusters %in% c(1,2,3))
# ggplot(data=as.data.frame(colData(sceset)),aes(x=CD135,y=CD150,col=cluster))+
# theme_bw()+
# geom_point()+
# scale_color_manual(values=c("grey80","red"))+
# #labs(color="stem+MEP+myeloid")+
# ggtitle("stem+MEP+myeloid")
ggplot(data=as.data.frame(colData(sceset))[sceset$cluster %in% c(4),],aes(x=CD150,y=CD135,col=cluster))+
theme_minimal()+
geom_point()+
xlim(3.5, 10)+
ylim(5, 13)+
scale_color_manual(values=darkcols[c(4)])+
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
ggtitle("stem")
ggplot(data=as.data.frame(colData(sceset))[sceset$cluster %in% c(2,4),],aes(x=CD150,y=CD135,col=cluster))+
theme_minimal()+
geom_point()+
xlim(3.5, 10)+
ylim(5, 13)+
scale_color_manual(values=darkcols[c(2,4)])+
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
ggtitle("stem+MEP")
ggplot(data=as.data.frame(colData(sceset))[sceset$cluster %in% c(2,3,4),],aes(x=CD150,y=CD135,col=cluster))+
theme_minimal()+
geom_point()+
xlim(3.5, 10)+
ylim(5, 13)+
scale_color_manual(values=darkcols[c(2,3,4)])+
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
ggtitle("stem+MEP+myeloid")
ggplot(data=as.data.frame(colData(sceset)[sceset$cluster %in% c(1,2,3,4),]),aes(x=CD150,y=CD135,col=cluster))+
theme_minimal()+
geom_point()+
xlim(3.5, 10)+
ylim(5, 13)+
scale_color_manual(values=darkcols)+
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
ggtitle("stem+MEP+myeloid+lymphoid")
dev.off()
RStudioGD
2
pdf("FACS_plot_per_clusters.pdf")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 1)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 1")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 1")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 1")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 2)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 2")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 2")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 2")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 3)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 3")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 3")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 3")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 4)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 4")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 4")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 4")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 5)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 5")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 5")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 5")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 6)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 6")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 6")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 6")
cluster = as.factor(pData(sceset)$sc3_7_clusters == 7)
ggplot(data=pData(sceset),aes(x=CD135,y=CD150,col=cluster))+geom_point()+ggtitle("cluster 7")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD127,col=cluster))+geom_point()+ggtitle("cluster 7")
ggplot(data=pData(sceset),aes(x=Sca1,y=CD16.32,col=cluster))+geom_point()+ggtitle("cluster 7")
dev.off()
y = DGEList(counts=expr_data)
y = estimateCommonDisp(y)
design <- model.matrix(~cell_state)
fit = glmFit(y,design=design)
lrt = glmLRT(fit,coef = 2)
cc_genes = rownames(topTags(lrt,n=Inf,p.value=0.05))
lrt = glmLRT(fit,coef = 3)
cc_genes = c(cc_genes, rownames(topTags(lrt,n=Inf,p.value=0.05)))
sce_hi = sceall[rownames(hvg.out),]
#sce_hi = sce_hi[!(rownames(sce_hi) %in% cc_genes),]
sce_hi <- sc3(sce_hi, ks = 4:10, biology = TRUE)
sce_hi <- sc3_estimate_k(sce_hi)
pData(sce_hi)$G1 = assignments$scores$G1
pData(sce_hi)$G2M = assignments$scores$G2M
sc3_plot_markers(sce_hi, k = 7,show_pdata = c(
"G1",
"G2M"
), p.val=0.01,auroc=0.75)
ggplotly(ggplot(data=NULL,aes(x=tsne_sc3dim1$Y[,1],y=rtsne_ph$Y[,1],col=as.factor(pData(sceset)$sc3_8_clusters)))+geom_point())
We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
ggplotly(ggplot(data=NULL,aes(x=tsne_sc3$Y[,1],y=tsne_sc3$Y[,2],col=rtsne_ph$Y[,1]))+geom_point())
We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
ggplotly(ggplot(data=NULL,aes(x=tsne_sc3$Y[,1],y=tsne_sc3$Y[,2],col=exprs(sceset)["Dach1",]))+geom_point())
We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
ggplotly(ggplot(data=NULL,aes(x=tsne_sc3$Y[,1],y=tsne_sc3$Y[,2],col=as.factor(pData(sceset)$sc3_8_clusters)))+geom_point())
We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
#plot_ly(x=tsne_sc3$Y[,1],y=tsne_sc3$Y[,2],z=rtsne_ph$Y[,1],color=as.factor(pData(sceset)$sc3_8_clusters))
plot_ly(x=tsne_sc3$Y[,1],y=rtsne_ph$Y[,1],color=as.factor(pData(sceset)$sc3_8_clusters))
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
plot_ly(x=tsne_sc3$Y[,2],y=rtsne_ph$Y[,1],color=as.factor(pData(sceset)$sc3_8_clusters))
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode